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Abstract 

A polemic arose recently about the applicability of the i-expansion method to the calculation of the ground 
state energy Eq of the Rabi model. For specific choices of the trial function and very large number of 
involved connected moments, the i-expansion results are rather poor and exhibit considerable oscillations. 
In this letter, we formulate the t-expansion method for trial functions containing two free parameters which 
capture two exactly solvable limits of the Rabi Hamiltonian. At each order of the t-series, E^ is assumed 
to be stationary with respect to the free parameters. A high accuracy of E^ estimates is achieved for small 
numbers (5 or 6) of involved connected moments, the relative error being smaller than 10^* (0.01%) within 
the whole parameter space of the Rabi Hamiltonian. A special symmetrization of the trial function enables 
us to calculate also the first excited energy Ei, with the relative error smaller than 10^^ (!%)■ 
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Introduction: This paper is about the i-expansion 
method of the calculation of low-lying energy spec- 
trum for quantum Hamiltonian systems. The ap- 
plication of the method to the Rabi model evoked 
some doubts about its reliability [ll l2| . In this let- 
ter, we propose such treatment of i-expansion se- 
ries which provides, in low approximation orders, 
extraordinarily accurate estimates of the ground 
state energy in the whole range of model's parame- 
ters. First we explain the t-expansion method, then 
summarize the variational approaches to the Rabi 
Hamiltonian, propose a stationarity treatment of 
the t-expansion with two-parameter trial functions 
and finally present accurate numerical results. 

t-expansion: The t-expansion technique, origi- 
nated by Horn and Weinstein [3], is a "series ex- 
tension" of the variational method. It is based on 
the following theorem. For any trial function |z/^), 
which has a non-zero overlap with the exact ground 
state of a Hamiltonian iJ, the function 



monotonously decays in t, approaching the 
ground-state energy Eq at asymptotically large t: 
limt^oo E{t) = Eq. The coefficients of the small-t 
expansion {/„} are known as the connected mo- 
ments. They can be expressed recursively in terms 
of the standard moments ^m — {ip\H™\ip) as 



J?r 



t^r, 



m-2 

E 

fe=0 



1 



k 



-^fc+l/^rn-fc-l 



(2) 



E{t)^ 






i-ty 



m=0 



im+l 



(1) 



for m > 2 and /i — fii. The t = estimate £'(0) = 
/i is the variational value, representing a rigorous 
upper bound for the ground state energy. 

Usually only a limited number of connected mo- 
ments can be evaluated. Since we are interested in 
the t — >■ oo behavior of E{t), one needs some ex- 
trapolation from the small-t series to large t. From 
among various schemes [J, la, Isl, |7| we choose the fol- 
lowing ones. The widely used Connected moments 
expansion (CMX) [J] considers E{t) to be a sum 
of exponentials. The estimate of Eq, available only 
for odd number of connected moments m = 2k + 1 
{k = 0,1,...), is given by 
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where the vector X^ = (/2, ..., Ik+i) and the matrix 
Tfc has elements (rfc)^ = h+j+i, {i,j = 1, • • ■ , k). 
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In particular, we have Eq = Ii, Eq ^ Ii— /I//3, 
etc. The second scheme, known in the hterature 
alSlillOl as the Canonical sequence method (CSM) 
6j, corresponds to a polynomial "deformation" of 
one exponential. The method is formulated in the 
inverse format, using the function t{E) instead of 
E{t); the series expansion of t{E) around E = Ii is 
deducible from ([T|) [5] . The estimate of the ground 
state energy, which involves m connected moments 
(to may be even or odd), reads 
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where a-^'t{Ii) means the fc-th derivative of t{E) 
at E = Ii. We have the same Eq and Eq as in 
CMX, E^*^ =h+ 2/1/3/(12/4 - ill), etc. 

Rabi Hamiltonian: The Rabi model 12| de- 
scribes the interaction between a bosonic mode with 
energy a; > and a two- level atom with the gap luq. 
Its Hamiltonian is 

H=]^ujoa'' +ujh'^b + g{a+ + a-){b'^ +b), (5) 

where g is the interaction constant, a^' , a are the 
Pauli matrices, 6^ and b are boson creation and an- 
nihilation operators, respectively. 

There exist two exactly solvable cases of the Rabi 
Hamiltonian. For g — the system decouples and 
the ground-state wavefunction is the tensor product 
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as the atom stays at the bottom level and the boson 
is in his lowest mode as well. The energy of this 
state is Eo = — wo/2. The other case is uo — 0, 
when the two atomic levels merge to a degenerate 
one [13|. The exact (two- fold degenerate) ground 
state 
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is the product of the eigenfunction of {a^ 
2ct^ and the coherent boson state 
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One ground state is specified by the parameters 
{x — 2g/uj,y — —1}, the conjugate one by the op- 
positely signed {x = —2g/uj,y = 1}. The ground 
state energy is Eo — —Ag'^/uj. 



The (trial) function ([7]) is intentionally written 
so that both x and y can serve as variational pa- 
rameters for an arbitrary ujq 7^ 0, a.s was done in 
the standard variational method [l3|, |lj| . The op- 
timized values of the parameters a:opt and yopt are 
determined by minimizing Eq {x,y) — Ii{x,y): 
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The interpolation between the two exact solutions, 
{x = 0,y = 0} at g = (O and the branch {x = 
2g/u!,y — --1} at luq = 0, is provided by the unique 
solution {a;opt,yopt} with components restricted to 
the intervals 



< Xopt < — , 



-1 < yopt < 0. 
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Bishop et al. [13| pointed out a conserved parity 
of the Rabi Hamiltonian, associated with the sign 
reversal transformation {x,y} — >■ {—x, ~y}. The 
trial function ([7]), which does not possess this sym- 
metry, will be referred to as non-symmetrized. Two 
symmetrized versions of ([7]) were proposed: 
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where the normalization constants are 
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(p) and (n) stand for the positive and negative par- 
ity, respectively. It can be shown that \i{j^P\x,y)) 
and \il)^"'^x,y)) are orthogonal to one another. 
Within the variational approach [13|, \'ip^P^{x,y)) 
implies the ground state energy Eq of parity (p) The 
application of the (n)-symmetrized trial function 
\ip^^>{x,y)) projects the ground state away and, 
consequently, implies the first excited energy Ei of 
parity (n). 

To calculate moments of the Rabi Hamiltonian 
with an arbitrary one of the three trial functions 
|^(a;, y)) or \tp'-P'^'{x, y)), we apply the commutator 
&6^ = I + b^'b and the useful formula 
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{k,l — 0,1,2,...) valid for each of four possi- 
bilities xi — ±x, X2 — ±2:. The moments 



with the non-symmetrized trial function /i,; 
{4>{x,y)\H"^\tp(x,y)) are found to be 
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etc. The moments with the (p) and (n) sym- 
metrized trial functions are obtained in the form 






8gxy 



ujg y'^ -1 

2 2/2 + 1 (l + y2)Vl-e-4-' 



a; a; 



2n±1 



(p.") 
M2 



- [2/2(coth; 

+ (tanhx2)±i] , 
i+^V-|-452(l + 2a;2) 
8gwa;y(l -t- 2x'^) 



(2/2 -}- 1)\/1 - e-4^' 



-^[y2(eothx2)±i 

-(tanhx2)±i] 
(8g2+w2)a;2 



[y2(cothj:2) 



2\±1 



-(tanhx ) 



2\±11 



(15) 



,(p.") 



etc. We calculated the moments fim and ^m up 
to m = 6. 

Motivation: The eigenstate of the Rabi Haniilto- 
nian with g = ^ was used as a trial function for 
the t-expansion in Ref. |l|. For the CMX extrap- 
olation with 5 connected moments, the numerical 
results for Eq are satisfactory only in the region of 
small g. Amore et al. [2] used the CMX scheme 
with up to 99 connected moments, without a real 
improvement of the previous results for intermedi- 
ate and large values of g; in some regions of model's 
parameters, they even encounter numerical instabil- 
ities (considerable oscillations) of the results. The 
authors conclude that the method is not reliable for 
practical purposes. 

The method: Our idea is to use the CMX 

and CSM versions of the i-expansion method 
with the Rabi variational trial functions, the non- 
symmetrized \ip{x,y)) ([7]) and the (p), (n) sym- 
metrized |'0*-^'"''(^; 2/)) dJ- Using the i-expansion 



.Er>{x,y), 



with m =^ 1,3,4,... connected moments involved, 
we have at disposal Eq which depends on free 
parameters {x,y}. In the lowest (variational) or- 
der m — 1, the optimized values of {x]^J^,y^p{} = 
{s^opt, 2/opt} are determined by the stationarity con- 
ditions © which imply the global energy minimum 
in the (x, y) space. The determination of the free 
parameters for to > 3 is based on the following ar- 
guments. Although the limit liiRm^oo Eq { 
if it exists, would not depend on x and y, our fi- 
nite truncations Eq do. If the free parameters are 
chosen properly in a convergence range of the series, 
Eq' converge smoothly to the exact i?o as to — ^ 00. 

In the opposite case, Eq oscillates quickly as to 
increases which is an indication of loss of conver- 
gence properties. To ensure at least a "local inde- 
pendence" of Eq on free parameters, we impose 
the stationarity conditions 11] 
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This equation determines {x ( ,y)j„^}. In contrast 

to the variational Eq , the optimized Eq is not 
a rigorous upper bound for the ground state en- 
ergy. Sometimes, there exist more solutions of Eqs. 
(J16p . It is obvious not to accept maxima and sad- 
dle points, but still we can have several minima 
and now the global one need not to be the best 
choice. In general, there exists only a unique curve 
(hypersurface) of optimized Xq™^ (wq, w,^), and of 

the coupled y^^{uJo,io,g), which is continuous in 
Rabi's parameters and simultaneously lies close to 
the variational solution; this will be our physical 
solution. All other non-physical solutions, forming 
disconnected "blind arms" , are ignored; details for 
specific cases will be given bellow. 

Numerical results: We apply both CMX method 
([3]) in TO = 1,3,5 orders and CSM method (HJ in 
TO, — 1,3,4,5,6 orders. The difference between 
the CMX and CSM results turns out to be very 
small. In overwhelming number of cases, the re- 
sults are slightly above the best ( "exact" ) estimates 
obtained by the straightforward diagonalization of 
the Hamiltonian matrix in an appropriate basis set 
|2l.ll3|. In the Rabi Hamiltonian, one parameter can 
be fixed as it merely sets the energy scale; we prefer 
to set Wo = 1. Then we choose some w and gradu- 
ally change g in the whole interval [0, 00]. Estimates 
of Eq {ujo = l,uj,g) are expected to be satisfactory 
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Figure 1: The plot of Xopt{g) for ojq ■ 
the non-symmetrized trial function (O. 
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for small and large values of g as the trial functions 
are, by construction, close to the exact solutions at 
g — and 5 — > cx) {luq = previously). The true 
problem is the region of intermediate values of g. 

First we present the results for the non- 
symmetrized trial function \ip{x,y)) ([7]). For any 
{uJo,uj,g}, the variational Eq possesses just one 
minimum at {xopt, 2/opt} lying in the interval ([T0|. 
Each of the functions a;opt(5) and yo-ptig) is contin- 
uous. The plot of Xopt(g) is represented in Fig. [1] 
by the solid curve [the picture is similar for yoptig)]- 
For 16(7^ < ijJqlo, the solution of Eq. (ITSl) giving the 

minimum of -Eq is trivial: {a;opt,J/opt} = {0,0}. 
This is why the curve Xo-pt{g) lies on the g axis up 
io g — 1/4. If g > 1/4, the solution becomes non- 
trivial and approaches the asymptotic Xopt = 2g/ijj 
(dashed line) for large g. Although this curve is 
continuous, it is non-analytic at the point g = 1/4. 
Going to Eq with tti > 3, the number of minimum 
solutions to Eq. (IT5|) can be larger than one in cer- 
tain intervals oi g. The minima curves can break at 
some points (beyond which there are no minimum 
solutions) or split into several curves. None of them 
goes continuously from small to large values of g. 

Nevertheless, the curves denoted by open dia- 
monds and circles in Fig. [T] can be used for small-g 
and large-5 cases, respectively, to obtain satisfac- 
tory results. For small 5, the trivial minimum {0, 0} 
can be directly inserted into all expressions. The 



Figure 2: The plot of Xop\,{g) for ujq = u) = 1, calculated 
with the (p) symmetrized trial function lllll l . The solid curve 
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the dashed line is the 



corresponds to the variational E^ 

large g-asymptotic 

parameters closest to the solid curve are denoted by ^, blind 

arms are denoted by A and X 



The CMX estimate e''^'^^: Optimized 



Table 1. Estimates of the energies Eo and Ei for 
the Rabi Hamiltonian with g = 5. 





UJQ = UJ = 1 


LJQ — l,UJ = 2 


p(6) 

/T^cxact 
^0 


-100.006250000 
-100.006265682 
-100.006265704 


-50.001250000 
-50.001262703 
-50.001262758 




-100.006250000 
-100.006265686 


-50.001250000 
-50.001262703 


variational Eg (0,0) — —luo/2 and 
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For g <t^ (jjQ^ this formula is consistent with the 
exact result with the relative error |^q (0,0) — 



ZT'cxact 
^0 



|/£^oxact Qf ^i-^^ Qj.jgj. ;^Q-5^ rpj^g large-5 re- 
sults, illustrated in the first window of Table 1, are 
even more precise. We present the variational result 
Eq , the CSM Eq and the numerically exact re- 
sult 0. The values of {a;opt,yopt} and {a;ipt,yopt} 
are very close to the asymptotic result {2g/uj, —1}. 
For example, for wq = cj — 1 and g — 5, the 
minimum of -B^*^^ is at {9.99997,-0.997364}. The 
formulas with symmetrized trial functions (see bel- 
low) lead to the results with comparable accuracy 
in both small-g and large-g regions. 

The results for the ground state energy obtained 
with the (p)-symnietrized trial function \ip^''''{x,y)) 
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Figure 3: The ground-state energy E^'^ for uiq = lj = g = 
1, calculated with the (p)-symmetrized trial function: The 
CSM results Q in m = I, 3,4, 5, 6 orders, the CMX results 
• in m = 1, 3, 5 orders. Dashed line is the numerically exact 
value [3. 



Figure 4: The CSM estimate of the ground state energy 
with the (p) symmetrized trial function, E^ '^ , for ujq = 1, 
various discrete values of ui and the interaction parameters 
g constrained to the interval [0, 1] . 



PT|) are presented in Fig. [5J The main advan- 
tage of the variational resuh Eq '^ in comparison 
with the non-symmetrized one is that the curve of 
unique minima, plotted as the solid curve, is not 
only continuous but also analytic for all g. In the 
considered higher orders Eq'^' (3 < m < 6) and 
for both CSM and CMX methods, there exists a 
unique counterpart of the variational curve of min- 
ima (represented by open diamonds) which is con- 
tinuous and free of singular points in the whole in- 
terval of g values; these are the accepted optimized 
minima. Similarly as for the non-symmetrized trial 
function, blind disconnected curves of minima ap- 
pear (open triangles and crosses); we ignore them. 
Since some of minima are very close to the varia- 
tional curve, the region of intermediate values of g 
is magnified in the inset of Fig. [5] 

First we discuss the results for small-g. Bishop 
et al. [13[ showed that for g <$^ luq the coordinates 



of minima for the variational E, 
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T5]) ] are, up to the term linear in g, given by 
2ff ., 2g 



•^opt 



a/w(wo + w) ' 



yopt 



,(i,p) 



Wo 



(18) 



Inserting these values into E^ and expanding in 

small g up to the g^ term we reproduce Eq. ([TT]). 

(3) 
derived from the non-symmetrized £^q . This co- 
incidence confirms that the non-symmetrized trial 
function gives adequate results not only for large g, 
where the corresponding formulas effectively merge 



because the corresponding x is large, but also in the 
region of small g. 

For Rabi's parameters luq = lo = g = 1, a quick 
convergence of the results for the ground state en- 
ergy Eq'^ to the exact value as m increases is 
shown in Fig. [3] The CMX data are represented 
by full circles, the CSM data by open circles; note 
that the CMX and CSM results are very close to 
each other. The (numerically) exact value [13| is 
represented by the dashed line. The improvement 
of the variational m — 1 result is remarkable al- 
ready for m — 3. As m increases, the convergence 
of the data to the exact value is excellent. 

Our procedure enabled us to calculate very 
quickly the CSM ground state energies Eq '^ for 
various sets of the Rabi Hamiltonian parameters, 
see Fig. 01 Without any loss of generality we set 
luq — 1. Each curve is labeled by the boson energy 
0.05 < w < 5. The interaction parameter g is con- 
strained to < g < 1. All values are correct within 
the resolution of the plots. The worst relative error 
of the order 10^ was achieved for medium values 
of 5. 

To document an extraordinary accuracy of the 
obtained results, in Fig. [5] we present the rel- 
ative error of the CSM estimates of the ground 
state energy E^ '^ [{p) symmetrized trial function, 
open circles] and of the first excited energy E\ '"'' 
[{n) symmetrized trial function, open triangles] for 
ujq = uji = 1 and an interval of g-values. Within 
the whole parameter space of the Rabi Hamilto- 




Figure 5; The relative error of the CSM estimates of the 
ground state energy Eq '^ [{p) symmetrized trial function, 

open circles] and the first excited energy _EJ '" [(n) sym- 
metrized trial function, open triangles], for luq = ui = 1 and 
an interval of g-values. 

Table 2. Estimates and relative errors of Ei for 
Wo = 1^ = 1 and g — 0.2. 

Estimate Rel. error 



El ' var. 


0.00324806 


0.39 


£;f '"^ CMX 


0.00233753 


0.000335 


£f '") CSM 


0.00234135 


0.00197 


TTicxact 


0.00233675 






nian, the relative error is smaller than 10 ^ for Eq 
and smaller than 10^^ for Ei. Like for example, 
our CSIVI results for the ground state energy Eq '^' 
are -0.69761396 for g = 0.3 and -0.87854267 
for g — 0.4, the (numerically) exact values Eq^'^'^^ 
are -0.69761529 for g = 0.3 and -0.87854932 for 
g = 0.4 [16|. The relative errors are 1.9 x 10^^ and 
7.8 X 10~^, respectively. As concerns the first ex- 
cited energy, Bishop et al. [13[ reported the largest 
relative error (almost 0.4) for the variational esti- 
mate E[ at the interaction constant g = 0.2. We 
see in Table 2 that this error goes down quickly in 
higher approximation orders. 

Conclusion: In conclusion, it turns out that the 
optimization of low orders of the i-expansion for the 
Rabi Hamiltonian improves remarkably the preci- 
sion of the variational ground state estimates. In 
the whole parameter range of Rabi model, with 
the symmetrized trial function the relative error is 
smaller than lO^'' (0.01%). The accuracy of the 
i-expansion method is not ensured if the trial func- 



tion is not properly chosen, as was seen in the case 
of non-symmetrized trial function and medium val- 
ues of the interaction constant g. 
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